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Ground-state magnetic properties of the diluted Heisenberg antiferromagnet on a square lattice are 
investigated by means of the quantum Monte Carlo method with the continuous-time loop algorithm. 
It is found that the critical concentration of magnetic sites is independent of the spin size S, and 
equal to the two-dimensional percolation threshold. However, the existence of quantum fluctuations 
makes the critical exponents deviate from those of the classical percolation transition. Furthermore, 
we found that the transition is not universal, i.e., the critical exponents significantly depend on S. 
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Since the discovery of high-temperature superconduc- 
tivity in cuprates §1, effects of non-magnetic impurities 
in two-dimensional Heisenberg antiferromagnets (HAF's) 
have been studied extensively. As is weW known, the an- 
tiferromagnetic (AF) long-range order in pure La2Cu04 
{S = 1/2 HAF) is immediately destroyed by a small 
amount of Sr substitution for La. Itinerant holes, doped 
into the Cu02 plane, cause magnetic frustration among 
Cu S = 1/2 spins and destroy the AF order. On the 
other hand, the AF long-range order is much robust 
against static magnetic impurities such as = Mg 
or Zn in La2Cu04, and Mg in K2C0F4 {S = 1/2 
Ising AF) or K2MnF4 (S=5/2 HAF). However, it is 
notable that the critical concentration of non-magnetic 
impurities strongly depends on the system in these 
cases. The critical concentration of La2Cui_a;Mga;04 and 
La2Cui_a:Znj;04 is observed Q as a; ^ 20%. It makes 
a sharp contrast with the others where the criti- 

cal concentration is almost equal to the two-dimensional 
percolation threshold {x ^ 40%). In La2Cui_2;Mg2;04 
and La2Cui_2;Znj,04, a quantum disordered phase may 
be realized due to strong quantum fluctuations. 

These experimental results, and also purely theoreti- 
cal interests, have motivated recent theoretical and also 
numerical studies j| [s] on the site-diluted HAF on the 
square lattice: 
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Here, — {Sf , Sf , S^) denotes the the quantum spin- 
s' operator at site i, and the nearest-neighbor coupling 
constant is antiferromagnetic (J > 0). The quenched di- 
lution factors {ci} independently take 1 or with proba- 
bility p and l—p, respectively, where p {= 1~ x) denotes 
the concentration of magnetic sites. 

In the classical case {S = 00), the present model at zero 
temperature is equivalent to the site-percolation prob- 
lem [S. It is well known that the system undergoes a 



second-order phase transition at the percolation thresh- 
old Pel, which is determined as 



Pel = 0.5927460(5) 



(2) 



Near pd, the stag- 
^ (p ~ Pci)'^"' with 




by the most recent simulation |10[| . 
gered magnetization vanishes as Ms 
/?^i = 5/36 = 0.13888- ■• §. 

The main subjects in the previous works have 
been focused on whether the critical concentration of the 
quantum systems (S < 00), referred to as p* hereafter, 
is identical to that of the classical model, pd, or not, 
and also on its dependence on the strength of quantum 
fluctuations specified by the spin size S. For S — 1/2, 
p* > Pel is suggested in these studies ^-||,|lll 
ample, p* = 0.655 and 0.695 are obtained in 
respectively. However, they are still inconclusive and crit- 
ical properties of the quantum phase transition and the 
quantum disordered phase (if exists) between p* and pd 
have not been well understood. 

In this paper, we report results of large-scale quantum 
Monte Carlo (QMC) simulations on the diluted HAF @) 
with 5* = 1/2, 1, 3/2, and 2. It is found, as we see below, 
that the AF long-range order at T = persists so long as 
a cluster of magnetic sites percolates, that is, p* = pd, 
even in the S = 1/2 case. However, we find the critical 
exponents of the quantum phase transition just at Pci are 
clearly different from those of the classical transition; the 
critical exponents vary depending on the value of S. 

In the present QMC simulation, use of the continuous- 
time loop algorithm [T^-[l5|] is crucial. It greatly reduces 
correlations between successive world-line configurations, 
and therefore makes it possible to perform highly reliable 
simulations on large lattices (up to L x L = 48 x 48) at 
extremely low temperatures (T ~ O.OOIJ). Another im- 
portant feature of the present algorithm is its ergodicity; 
the winding number of world lines around vacant sites 
can change and the ground-canonical ensemble can be 
simulated. We use the periodic boundary conditions. For 
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FIG. 1. System-size dependence of Sb{L,T = 0,p)/L'' in 
the case of S = 1/2 at p = 1, 0.875, 0.8125, 0.75, 0.6875, 
0.65, and 0.625 (from upper to lower). The dashed lines are 
obtained by least-squares fitting for the largest three system 
sizes for each p. The extrapolated values are denoted by solid 
symbols. The data with p < 0.6875 and L > 24 are also 
shown in the inset. 



each sample, 10^ Monte Carlo steps (MCS) are spent for 
measurement after 10^ MCS for thermalization. At each 
parameter set {L,T,p), physical quantities are averaged 
over 100 - 1000 samples. 

First, we discuss the zero-temperature staggered mag- 
netization Ms{p) for p > Pel, which is calculated as 



lim lim 



35s(£,r,p) 



in terms of the static structure factor defined by 
1 - 



Ss{L,T,p) 



(3) 



(4) 



at the momentum k = (tt, tt) , where d is the spatial di- 
mension {d — 2). The bracket in Eq. denotes both of 
the thermal average and the average over samples. In the 
present simulation, we use an improved estimator to cal- 
culate Ss{L, T,p), by which the variance of data is greatly 
reduced. 

At temperatures lower than the gap of a finite system, 
Ss{L,T,p) converges to its zero-temperature value quite 
rapidly (probably exponentially). To obtain Ss{L,T — 
0,p) at each p and L, we perform QMC simulations at 
low enough temperatures so that Ss{L,T,p) exhibits no 
temperature dependence besides statistical errors. Note 
that as p decreases, the gap of a system of linear size 
L becomes smaller, and therefore lower temperature is 
needed For the S ^ 1/2 case, T is taken as 0.002 J 
at p = 0.625 and L = 48. 

In Fig. g, we plot S',(L, 0,p)/L'* against 1/L for 0.625 < 
p < 1 in the 5* = 1/2 case. It is clearly seen that the data 
at each concentration fall on a straight line for large L. 
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FIG. 2. Concentration dependence of the (normalized) 
staggered magnetization in the cases of S = 1/2 (square) and 
1 (triangle). The dotted lines are guides to eyes. In the inset, 
we show the double-logarithmic plot of the staggered magne- 
tization against (p — pd)- The dashed lines are obtained by 
least-squares fitting for p < 0.70. For comparison, the stag- 
gered magnetization in the classical limit is also plotted by 
the solid line and diamonds. 



In the clean system {p = 1), the leading finite-size correc- 
tion is shown to be of 0{1/L) according to the spin- wave 
theory jl^] . The similar behavior for p < 1 indicates that 
there exists an AF long-range order with massless exci- 
tations, such as spin waves, even in the presence of im- 
purities. Thus, the staggered magnetization in the ther- 
modynamic limit is obtained by linear extrapolation in 
1 /L for three largest system sizes at each p. 

The final results for the zero-temperature staggered 
magnetization are shown in Fig. |^. It is seen that the 
staggered magnetization remains finite even at p = 0.625 
both in the = 1/2 and 1 cases. The possibility that p* is 
greater than 0.625 is excluded definitely, and the behavior 
of the whole magnetization curve strongly suggests that 
the critical concentration is identical to the percolation 
threshold {p* = pc\). Actually, the staggered magneti- 
zation seems to vanish algebraically towards p = pd as 
seen in the inset of Fig g. By least-squares fitting, we 
estimate the critical exponent f3 as 0.46(3) and 0.32(3) 
for S = 1/2 and 1, respectively. Our conjecture, p* = pci, 
is also supported by the power-law behavior of the zero- 
temperature static structure factor just at p = pd'. 



(5) 



with * = 1.17(6) and 1.57(3) for S = 1/2 and 1, respec- 
tively (Fig. |). 

The critical exponents f3 and "if, obtained by the above 
analysis, differ from those of the percolation model {/3 — 
5/36 and ^E" = 43/24 |^), that is, the universality class of 
the quantum phase transition at p = pd is different from 
the classical one. Furthermore, we found that the value of 
the critical exponents depends on S. This means that not 
only the fractal nature of the lattice geometry, but also 
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FIG. 3. System-size dependence of Ss{L,T = 0,p) at 
p — Pel in the cases of S — 1/2 (square) and 1 (triangle). The 
dashed lines are obtained by least-squares fitting for L > 24. 
The slope of the lines gives the exponent <& (Eq. (^)). The 
data in the classical limit are also plotted by diamonds. 



the existence of quantum fluctuations and their strength, 
or the value of S, are relevant to the criticality. Hence, 
it is natural to introduce the two following assumptions 
on the spin correlation on spin clusters. 

At p = Pel, all clusters are fractal with a fractal di- 
mension D (in two dimensions, D = 91/48) First, we 
assume that the staggered spin correlation between two 
sites, i and j, on a cluster behaves as 



C{i,j)--rJ forr,j»l, 



(6) 



where 



Here, we introduce a new S- 



dependent exponent a ~ a(S) > 0. In the classical case, 
C(i,j) takes a constant independent ofrij, i.e., 0(00) — 
0. Together with the cluster-size distribution dA, p — pd, 
predicted by the percolation theory [p|, we obtain 



2D-d-a 



(7) 



where the first summation on the first line is taken over 
clusters on the lattice. Comparing Eq. (0) with Eq. (||), 
we obtain a scaling relation: 

TABLE I. Summary of critical exponents 'i, z, 0, and u. 

and z are obtained by the ESS shown in Fig. |4 and P 
is estimated from the staggered magnetization (FigTg). The 
exponent 1/ is calculated from (3 and ^ thorough the scaling 
relations, Eqs. (||) and (pT|). 
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1.27(2) 


2.54(8) 
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1.2(1) 
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1.57(3) 


1.58(10) 


0.32(3) 


1.5(2) 
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1.60(3) 


1.55(10) 
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1.69(7) 


1.31(20) 
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1.79167 




0.13889 


1.33333 



On the other hand, for p > pci, the percolation the- 
ory 1^ says that there exists a characteristic length A(p), 
below which the percolated cluster exhibits a fractal na- 
ture similar to that of clusters sX p — pd, but it has 
the ordinary dimension {d = 2) at longer length scale. 
For L ^ A(p), Ss{L, T,p) is dominated by the percolated 
cluster. We introduce the second assumption that there 
exist no other macroscopic characteristic lengths except 
X{p) even in the S < 00 cases. Since there exists an AF 
long-range order as shown before, the correlation func- 
tion is expected to obey the following scaling form pM : 



C{i,j;p) 



A(p)-° 



for nj » A(p) , 
— » 0, so that Eq. ( 



(9) 



where C{x) = const, for a; — » U, so that Eq. (|Bp is repro- 
duced for 1 <C rij <^ A(p). 

Using the fact that the characteristic length scale X{p) 
diverges as {p — Pd)^" near pd with v — 4/3 and as- 
suming the ordinary algebraic temperature dependence, 
one finally reaches a full finite-size scaling (FSS) form of 
the structure factor: 



Ss{L,T,p)^L 



2D-d 



-"Ss{L'^''{p~Pd),L'T), (10) 



where the ^'-dependent scaling function Ss{x,y) takes a 
constant at {x,y) = (0,0), and 5s(a;,0) ~ a;-(2£'-2d-a),. 
for a; 3> 1. In Eq. (|lO|), we introduce the dynamical ex- 
ponent z, which relates the energy scale with the length 
scale @. In the present case, 2; > 1 is expected, since the 
Lorentz invariance is broken due to the existence of im- 
purities. Although the algebraic T dependence assumed 
in Eq. ( |l0|) is not guaranteed a priori, finite-temperature 
data of the staggered structure factor at p = pd is well 
scaled for S = 1/2, 1, 3/2, and 2, by using our FSS form 
(l^) with the exponents 4" and z listed in TABLE |, as 
shown in Fig. ^. The values of ^' are consistent with those 
obtained by the FSS at T = (Eq. ^ and Fig. |). Note 
that not only "if but also z depend on S. In addition, "if 
seems to approach the classical value monotonically as S 
increases. 

The critical exponent /3 is also expressed in terms of 
a. By taking L — > 00 after taking the limit of T ^ in 
Eq. (p^), one obtains 



2/3 = -{2D - 2d - a)iy . 



(11) 



The values of v calculated from \1/ and f3 through Eqs. (||) 
and ((ll| ) are satisfactorily consistent with ^ = 4/3 (TA- 
BLE 1)7 which supports the validity of the scaling argu- 
ment discussed above. 

In summary, we have investigated the ground-state 
phase transition of the diluted HAF with S = 1/2, 1, 
3/2, and 2. Contrary to the previous works 
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FIG. 4. Scaling plot of 5's(i, T,p) at p = pd for (a) S = 1/2 
and 1, and for (b) S = 3/2 and 2. 



our present QMC study has shown that the critical con- 
centration is equal to the classical percolation threshold 
even in the S = 1/2 case. Concerning the relation to 
experiments, it becomes clear that the present model (|l|) 
is too much idealized to predict the magnetic proper- 
ties observed in the real materials. Other effects, such 
as next-nearest-neighbor interaction, should be included 
properly into the model Hamiltonian. On the other hand, 
in the theoretical point of view, it has been revealed that 
the present model contains quite rich physics; the critical 
exponents vary depending on S. We have introduced a 
scaling form with two non-classical exponents, a and z, 
which consistently explains the behavior of the simula- 
tion data together with the exponents D and v of the 
classical percolation transition. To our best knowledge, 
this is the first time that such an S'-dependent quantum 
critical behavior is discovered. 

Most of numerical calculations for the present work 
have been performed on the CP-PACS at University 
of Tsukuba, Hitachi SR-2201 at Supercomputer Center, 
University of Tokyo, and on the RANDOM at Materi- 
als Design and Characterization Laboratory, Institute for 
Solid State Physics, University of Tokyo. The present 
work is supported by the "Large-scale Numerical Sim- 
ulation Program" of Center for Computational Physics, 
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the Future Program (Computational Science and Engi- 
neering)" of Japan Society for the Promotion of Science. 
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Education, Science, Sport and Culture of Japan. 
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